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Abstract 

Goal-oriented error estimates (GOEE) have become popular tools to quan- 
tify and control the local error in quantities of interest (Qol), which are often 
more pertinent than local errors in energy for design purposes (e.g. the mean 
stress or mean displacement in a particular area, the stress intensity factor 
for fracture problems,...). These GOEE are one of the key unsolved problems 
of advanced engineering firms in the aerospace industry. Residual-based error 
estimators can be used to estimate errors in quantities of interest for finite 
element approximations. This work presents a recovery-based error estima- 
tion technique for Qols whose main characteristic is the use of an enhanced 
version of the Superconvergent Patch Recovery (SPR) technique previously 
used for error estimation in the energy norm. This enhanced SPR technique is 
used to recover both the primal and dual solutions. It provides a nearly stat- 
ically admissible stress field that results in accurate estimations of the local 
contributions to the discretization error in the Qol and, therefore, in an ac- 
curate estimation of this magnitude. This approach leads to a technique with 
reasonable computational cost that could be easily implemented into already 
available finite element codes. 

KEY WORDS: goal-oriented, error estimation, recovery, quantities of interest, error con- 
trol, mesh adaptivity 

1 Introduction 

Verification and quality assessment of numerical simulations is a critical area of 
research. Techniques to control the error in numerical simulations and verify the 
approximate solutions have been a subject of concern for many years. Further 
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development is expected to have a profound impact on the reliability and utility of 
simulation methods. Effective methods for assessing the global discretization error 
in energy have been extensively developed since the late 70s. However, controlling 
the error in quantities of interest to engineers, e.g. the average stress or average 
displacement in a given region or the stress intensity factors is not equivalent to 
controlling the error in energy. Yet, the widely stated need of practicing engineers 
is to evaluate and minimise the discretisation error in these "quantities of interest." 
Less attention, however, has been given to the emerging, more advanced, "goal- 
oriented" error estimates. This paper is a step in this direction which provides a 
simple procedure to evaluate engineeringly-relevant quantities of interest. 

Energy-based error estimates can be classified into different families [1,2]: residual- 
based error estimators, recovery-based error estimators, dual techniques, etc. Residual- 
based error estimators originally introduced by Babuska and Rheinboldt [3] are char- 
acterised by a strong mathematical basis. To estimate the error, they consider local 
residuals of the numerical solution. By investigating the residuals occurring in a 
patch of elements or even in a single element it is possible to estimate the errors 
which arise locally. These methods were improved with the introduction of the resid- 
ual equilibration by Ainsworth and Oden [1]. Recovery-based error estimates were 
first introduced by Zienkiewicz and Zhu [4] and are often preferred by practitioners 
because they are robust, simple to use and provide an enhanced solution which 
may be used as output. Further improvements were made with the introduction 
of new recovery procedures such as the superconvergent patch recovery (SPR) pro- 
posed by the same authors [5, 6] and many papers following [7, 8, 9, 10, 11, 12]. 
Recovery techniques were extended to enriched approximations in [13, 14, 15, 16] 
and to Hellinger-Reissner smoothing-based finite elements in [17] and the role of 
enrichment and statical admissibility in such recovery procedures discussed in [18]. 
Duality-based techniques rely on the evaluation of two different fields, one compat- 
ible in displacements and another equilibrated in stresses [19]. 

Most techniques used to obtain error estimates prior to the mid 90s were aimed 
to evaluate the global error in the energy norm. Nevertheless, as discussed above, 
the goal of numerical simulations is generally not the determination of the strain 
energy alone, but the reliable prediction of a particular quantity of interest (Qol) 
which is needed for making decisions during the design process. Therefore, it is 
important to guarantee the quality of such analyses by controlling the error of the 
approximation in terms of the Qol rather than the global energy norm. Significant 
advances in the late 90s introduced a new approach focused on the evaluation of 
error estimates of local quantities [1, 20, 21]. In order to obtain an error estimate for 
the Qol two different problems are solved: the primal problem which is the problem 
under consideration, and the dual or adjoint problem which is related to the Qol. 
The adjoint problem is the same as the primal problem except for the loads and 
boundary conditions, which are determined as a function of the Qol. Goal-oriented 
error estimators have been usually developed from the basis of residual formulations 
and the widely used strategy of solving the dual problem [1, 22]. Although limited, 
the use of recovery techniques to evaluate the error in quantities of interest can be 
found in [23, 24], and considering dual analysis in [25]. 

In this paper we propose to extend the recovery techniques presented in [12, 26, 
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15, 27] to evaluate accurate error estimates for linear Qols. The recovery technique 
requires the analytical expressions of the body loads and boundary tractions which 
must be evaluated also for the dual problem in order to consider equilibrium con- 
ditions. In Section 2 we define the problem and the finite element approximation 
used. Section 3 focuses on the representation of the error in the Qols and the solu- 
tion of the dual problem. Section 4 deals with the definitions of the dual problem 
and describes the expressions of body loads and boundary tractions required for the 
stress recovery of the dual problem for different linear Qol. In Section 5 we introduce 
the nearly equilibrated recovery technique used to obtain enhanced stress fields for 
both the primal and dual solutions. In Section 6 we present some numerical results 
and, finally, conclusions are drawn in Section 7. 

2 Problem statement and solution 

In this section, we briefly present the model for the 2D linear elasticity problem. 
Denote, in vectorial form, cr and e as the stresses and strains, D as the elasticity 
matrix of the constitutive relation cr = Ds, and the unknown displacement field u, 
which values in f2 C R 2 . u is the solution of the boundary value problem given by 

-Vcr(u) = b infi (1) 

a(u)-n = t on IV (2) 

u = u on r c , (3) 

where Tjy and Tjj denote the Neumann and Dirichlet boundaries with dQ = 
r^UT^i and IV fl = 0, b are body loads and t are the tractions imposed 
along Ttv- 

Consider the initial stresses cr and strains Eq, the symmetric bilinear form a : 
(V + u) x V — > R and the continuous linear form £ : V — > R are defined by 

a(u,v):= f <r T (u)£(v)dfi = f cr T (u)D- 1 cr(v)dfi (4) 
Jn Jn 

e(v) := [ v T bdfi + I v T tdr + / <r T (v)£ dfi - / £ T (v)<r da (5) 
Jn ir N Jn Jn 

With these notations, the variational form of the problem reads [28] 

Find u G (V + u) : Vv G V a(u, v) = £(v) (6) 

where V is the standard test space for the elasticity problem such that V = {v | v G 

[#W,v|r D (x) = 0}. 

Let u h be a finite element approximation of u. The solution for the discrete 

counterpart of the variational problem in (6) lies in a subspace (V h + u) C (V + u) 

associated with a mesh of finite elements of characteristic size h, and it is such that 

Vv ft G V h C V a(u h , v h ) = l(y h ) (7) 
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3 Error control 

3.1 Error representation 

In this section we define a general framework for goal-oriented error analysis. First, 
we assume that the finite element discretization error is given by e = u — u h . To 
quantify the quality of u h in terms of e different methods have been proposed, 
generally based on the energy norm ||e|| induced by a/ a(e, e) and written as 



where cr and cr represent the exact and finite element stress fields, respectively. 
Following Zienkiewicz-Zhu [4], an estimate of the error ||e es || can be formulated in 
the context of FE elasticity problems by introducing the approximation 



where cr* represents the recovered stress field, which is a better approximation to 
the exact solution than cr h . Similarly, local element contributions can be evaluated 
from (9) considering the domain Q e of element e. 

The error estimate measured in the energy norm as given in (9) can overestimate 
or underestimate the exact error although it is asymptotically exact 1 provided that 
cr* has a higher convergence rate than the FE solution. In order to improve the 
quality of our recovery-based estimate it is useful to include the error bounding 
property. Recently, Diez et al. [26] and Rodenas et al. [27] have presented a 
methodology to obtain practical upper bounds of the error in the energy norm 
||e|| using an SPR-based approach where equilibrium was locally imposed on each 
patch. The recovered stresses obtained with this technique were proven to provide 
very accurate estimations of the error in the energy norm. 

Note that we have followed here the definition of error estimates provided in 
References [4, 7, 29] and that we reserve the term error bounding for techniques which 
provide upper or lower bounds of the error. The semantics used in Mathematics 
refer to (9) as an error indicator and any technique providing bounds as an error 
estimator. 

3.2 Error in quantities of interest: duality technique 

In this section we show how error estimators measured in the energy norm may be 
utilised to estimate the error in a particular quantity of interest [1]. The strategy 
consists in solving a primal problem, which is the problem at hand, and a dual 
problem useful to extract information on the quantity of interest (identical to the 
primal problem except for the applied boundary conditions and internal loads). 
From the FE solutions of both problems it is possible to estimate the contribution 
of each of the elements to the error in the Qol. This error measure allows to adapt 

1 the approximate error tends towards the exact error 
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the mesh using procedures similar to traditional techniques based on the estimated 
error in the energy norm. 

Consider the linear elasticity problem given in (6) and its approximate FE solu- 
tion u h . For the sake of simplicity let us assume homogeneous Dirichlet boundary 
conditions in this problem. This problem is related to the original problem to be 
solved, that henceforth will be called the primal problem. 

Let us define Q : V — > R as a bounded linear functional representing some 
quantity of interest, acting on the space V of admissible functions for the problem 
at hand. The goal is to estimate the error in functional Q when calculated using 
the value of the approximate solution u h as opposed to the exact solution u: 

Q(u)-Q(u h ) = Q(u-u h ) = Q(e). (10) 

As will be shown later, Q(v) may be interpreted as the work associated with a 
displacement field v and a distribution of loads specific to each type of quantity of 
interest. If we particularize Q(v) for v = u, this force distribution will allow to 
extract information concerning the quantity of interest associated with the solution 
of the problem in (6). 

A standard procedure to evaluate Q(e) consists in solving the auxiliary dual 
problem (also called adjoint or extraction problem) defined as: 

Find w Q G (V + w Q ) : Vv G V a(v, w Q ) = Q(y). (11) 

An exact representation for the error Q(e) in terms of the solution of the dual 
problem can be simply obtained by substituting v = e in (11) and remarking that 
in case wq = 0, for all Wq G V h , due to the Galerkin orthogonality, a(e, Wq) = 
so that 

Q(e) = a(e, w Q ) = a(e, w Q ) - a(e, w£) = a(e, w Q - w£) = a(e, e Q ). (12) 
Therefore, the error in evaluating Q(u) using u h is given by 

Q(u) - Q(u h ) = Q(e) = a(e, e ) =[(*„- <rj) D" 1 (a d - <r*) dfi (13) 

Jn 

where a v is the stress field associated with the primal solution and cr^ is the one 
associated with the dual solution. 

Following [24] we can introduce a first coarse upper bound for the global error 
in the Qol employing the Cauchy-Schwarz inequality that results in: 

Q(e) = a(e,e Q )<||e||||e Q || (14) 

Then, the evaluation of the error in the Qol is now expressed in terms of the error 
in the energy norm for the dual and primal solutions, for which several techniques 
are already available. However, this upper bound is rather conservative due to the 
use of the Cauchy-Schwarz inequality [24]. Moreover, it does not provide a local 
indicator that can be used to guide the adaptivity process. 
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For a discretization with n e elements let us consider the local test space V e = 
{v|n e : v G V} and a e : V e x V e — > R as the associated bilinear form a to an element 
Q e : Q = |J We Q e such that 

Vu, vGV a(u, v) = a e (u, v) (15) 

n e 

To obtain sharper error measures and local indicators, we can decouple the ele- 
ment contributions as shown in [24] such that: 

(5(e) = a(e,eg) = ^a e (e,e Q ) from(12) (16a) 
<^|a e (e,e )| (16b) 
<5>|k||e Q |k (16c) 

n e 

< l|e||||e Q || (16d) 

Note that in (13) and the indicators derived afterwards, the error in the Qol is 
related to the errors in the FE approximations u h and Wq. On that account, any of 
the available procedures to estimate the error in the energy norm may be considered 
to obtain estimates of the error in the quantity of interest. 

Similarly from (8) to (9), we can derive an estimate of the error in the Qol from 
(13) which reads 

Q(e) « Q(e es ) = f (*• - <rj) D 1 (<r* - a h d ) dfi (17) 
Jn 

where cr* and a* d represent the recovered stress fields for the primal and dual prob- 
lems. Here, we expect that the more accurate the recovered stress fields cr* and cr*,, 
the more accurate the estimation of the local contributions to the global error (this 
is of special interest in goal-oriented adaptivity) and the sharper the error estimate 
in the quantity of interest. 

In order to obtain accurate representations of the exact stress fields both for the 
primal and dual solutions, we propose the use of the locally equilibrated recovery 
techniques introduced in [12, 15, 27] and described in Section 5. This technique, 
which is an enhancement of the superconvergent patch recovery (SPR) proposed in 
[5], enforces the fulfilment of the internal and boundary equilibrium equations locally 
on patches. For problems with singularities the stress field is also decomposed into 
two parts: smooth and singular, which are separately recovered. 

Continuing with the ideas in [24], and using the recovery technique to obtain 
accurate estimates of the error with the decoupled approach, we can formulate using 
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(16) the following error indicators: 



<2(e) w Q(e es ) = Si = ^a e (e es ,e Qies ) 



(18a) 



<S 2 : 




(18b) 



Yl\\ e ™\\n e \\ e Q, 



,es||n e 



(18c) 



< S 4 : 



e, 



(18d) 



The properties of these error estimates are discussed in [24] in more detail. In a 
general sense, the idea is that the more accurate the estimate, the less guaranteed is 
the upper bound property. The only guaranteed upper error bound is given by (18d), 
under the assumption that the recovery technique correctly enforces equilibrium 
conditions and provides upper error bounds in the energy norm for both the primal 
and dual solutions [26]. These indicators are shown in the numerical results for 
comparison. 

4 Analytical definitions for the dual problem 

The recovery procedure based on the SPR technique denoted as SPR-CX, and de- 
scribed in [12, 15, 26, 27], relies on the a priori knowledge of functions b and t to 
impose the internal and boundary equilibrium equations. For the primal problem 
these values (b p and t p ) are already at hand, as seen in (1, 2), but for the dual 
problem they are not readily available. In this Section we present a compilation of 
procedures to obtain the analytical dual loads for a few typical linear quantities of 
interest, such as the mean values of displacements and stresses in a sub-domain of 
interest Qj. Basic compilations have been presented in [30, 31] and, in the same line, 
in [32], where the authors defined relations between the natural quantities of interest 
and dual loading data. Here, we have also included other cases like the evaluation 
of dual loads when the quantity of interest is the generalized stress intensity factor 
(GSIF). 

4.1 Mean displacement in Qi 

Let us assume that the objective is to evaluate the mean value of the displacements 
along the direction a in a sub-domain of interest Qi C Q. The functional for the 
quantity of interest can be written as: 



where \Qi\ is the volume of Qi and c Ua is a vector used to select the appropriate 
combination of components of u. For example, c Ua = {1,0} T if a is parallel to the 
x-axis. Now, considering v G V in (19) results in: 




(19) 




(20) 
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Note that the term c Ua /|fij| formally corresponds to a vector of body forces in 
the problem defined in (6). Therefore, we can interpret = c Ua /|fij| as a constant 
vector of body loads applied within the sub-domain of interest fli which can be used 
in the dual problem to extract the mean displacements. 



4.2 Mean displacement along I\ 

For the case where the quantity of interest is the functional that evaluates the mean 
value of the displacements along a subset Tj of the Neumann boundary the 
expression reads: 

Q(u) = u a \ Vi = / u T c Ua dT (21) 

being the length of Tj and c Ua a Boolean vector used to select the appropriate 
component of u. Again, considering v e V in (21) we have: 

Note that the term c Ua /|rj| can be interpreted as a vector of tractions applied 
along the boundary in the problem defined in (6). Thus, = c Ua /|rj| is a vector 
of tractions applied on Tj and that can be used in the dual problem to extract the 
mean displacements along IV 



4.3 Mean strain in Q, 



In this case we are dealing with the first derivatives of displacements. Consider for 
instance the component e xx . This term can be expressed as: 



du x 
dx 



d_ d_ 
dx dy 



1 




V-(C xx u) (23) 



Operating similarly for the other components of the strain vector we would have: 
e yy = V ■ (C ra u) e xy = V • (C xy u) (24) 
Taking (24), the mean value of the strain e xx in Qj can be evaluated as: 



Q(u) = e xx \ ni 



1 



i a 



1 



v • (c xx u) an 



(25) 



Hi 



Assuming that u is such that we can apply the divergence theorem the previous 
equation can be written as: 



Q(u) 



|n 4 = J (C xx u) T ndT 



(26) 



where n is the unit outward normal vector. The domain integral is then trans- 
formed into a contour integral. Finally, choosing v e V instead of u in the previous 
expression and reordering we obtain: 



dr 



v T t d dr 



(27) 
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Note that the term (C^n) / formally corresponds to a vector of boundary 
tractions along the boundary of Therefore, we can interpret = (C^n) / 
as a vector of boundary tractions applied along the boundary of the sub-domain of 
interest Tj, which can be used to extract the mean strain in 

For any given linear combination of strains e v , denoted as 

e v = c €xx e xx + c tyy e yy + c txy e xy = cje = V • (Gu) (28) 

where 

G = C tj , x C XX + C eyy Cyy + C exy C X y (29) 

we could evaluate an expression similar to (27) such that 

0(v) = jf v T (j^j dr = jf v T t d dr. (so) 

4.4 Mean stress value in 

Let us consider now as Qol the mean stress value a a |n 4 given by a combination of 
the stress components in a domain of interest which reads: 

Q(U) = ffalfli = 7^7 / ^adfi = 7^7 I C^Crdtt. (31) 
</Qi |"i| 

The stress components are linear combinations of the strain components. Thus, 
any combination of stress components a a can be expressed as: 

a a = c T c * = <De = c£e, (32) 

where c Ua is the vector used to define the combination of stress components and c e 
defines its corresponding combination of strain components. Therefore, one could 
evaluate the mean value of any combination c CTq of stress components in f2j using 
(30), where G is evaluated using the coefficients of = in (29). 

4.5 Mean stresses and strains in a domain of interest from 
the perspective of initial stresses and strains 

Previously, the expressions to calculate mean stresses or strains in a domain of 
interest Qi have been obtained using the divergence theorem. However, we can 
calculate these quantities of interest considering them as an initial stress load case 
in the dual problem. This can be proved immediately 

Q(u) = <r a | ni = t^t / <£>dn = / §r\<rdn (33) 

I 1 1 J ^li J fij I 1 1 

From (5) we can define e o d = /|fij| corresponding to the term of initial 
strains that we need to apply in the dual problem to extract the value of a a \n r A 
similar formulation can be derived for the case of the mean strains in Qi such that 

<*0,d = C IJ\ Q i\- 
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4.6 Mean tractions along I\ 

Let t = {t n ,t t } T , with t n and t t the normal and tangential components of the 
tractions vector t. Let us assume that we want to evaluate, for example, the mean 
normal tractions along boundary Tj. The functional that defines the mean tractions 
along the boundary Tj can be expressed as 

Q(u) = i n = / t T cdr (34) 

Using (34) and considering the extraction vector c and the rotation matrix Rp 
that aligns the tractions normal to the boundary Tj we have: 




= / {t x t y }u d dT (35) 

In (35) the term = Rpc/|Fj| corresponds to a vector of displacements used as 
Dirichlet boundary conditions for the dual problem used to extract the mean value 
of the normal tractions along IY 

Note that in this case Wq = u d ^ 0, then (12) does not hold. We redefine the 
dual problem in (11), Vv e V, such that: 

a(v, Wo) =0 inO 
w Q = u d on Ti 



The dual solution can be expressed as wq = Wg + wq, where w Q = on Tj. 
Assuming that wq = u d is in the FE solution space, the FE approximation for (36) 
is also decomposed into two parts Wg = Wq° + wq where, again, Wq° = on Tj. 
Therefore, we have for the dual problem 

Vv G V a(v,w Q ) = a(v,w° Q ) + a(v,w Q ) = 
VvGF a(v, Wq) = -a(v, w q ) 

Substituting v = e in (37), using the Galerkin orthogonality property, a(e, Wq°) = 
0, and considering that eg = e Q we write: 

«( e > w o - w o°) = «( e > e Q) = «( e > e o) = -«( e > w q) (38) 
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Similarly, the Qol can also be rewritten by means of the divergence theorem. 
Thus, generalizing (35) Vv G V we have: 

Q(v) = / (Q<r(v)) r u d dr = / <r(v) T e(u d )dfi = a(v,u d ) (39) 

where Q is the matrix form of the Cauchy's law considering the unit normal vector 
to Tj, such that t r . = Qcr. Thus, Q(e) = a(e, u d ) = a(e, wq) and substituting in 
(38) we obtain the error for this Qol: Q(e) = — a(e, e Q ). 



4.7 Generalised stress intensity factor in 



Consider the problem of evaluating the generalised stress intensity factor (GSIF), 
that characterises the singular solution in problems with reentrant corners and 
cracks, as the quantity of interest. 

From [33, 34] we take the expression to evaluate the GSIF, which reads: 



dq 

dx. 



(40) 

'rii v / ux j 

where are the displacement and stress fields from the FE solution, u^ 2 \ a^> 

are the auxiliary fields used to extract the GSIFs in mode I or mode II and C is a 
constant that is dependent on the geometry and the loading mode, q is an arbitrary 
C° function used to define the extraction zone which must take the value of 1 at the 
singular point and at the boundaries and Xj refers to the local coordinate system 
defined at the singularity. 

Rearranging terms in the integral in (40) we obtain: 



Q(u) = K (1 ' 2) 



C 



(2) 

u\ q,i 

(2) 

u\ q o 

u (2) a \ u (2) a 
u 2 g 5 i + v,! q,2 



(1)\T 



1 

C 



(2) (2) 

°11 9,1 + °21 9,2 

(2) (2) 

°12 9,1 + a 22 9,2 



dtt (41) 



which could be rewritten as a function of initial strains e and body loads b: 

Q(u) = = [ (cr( 1 )(u)) T e + (u^fbdQ (42) 

Thus, if we replace u with the vector of arbitrary displacements v, the quantity 
of interest can be evaluated from 

Q(v) = / Lv T De 0i(i dfi + f v T b d dfi. (43) 

Hence, the initial strains and the body loads per unit volume that must be 
applied in the dual problem to extract the GSIF are defined as 



£o,d — 



1 

C 



(2) 

«i 9,1 

(2) 
«2 9,2 
(2) , (2) 



C 



(2) (2) 
^11 9,1 + ^21 9,2 

(2) (2) 
°12 9,1 + ^22 9,2 



(44) 
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5 Recovery technique 



As noted in previous sections, a widely used technique to control the error in the 
energy norm in the finite element discretization is the Zienkiewicz-Zhu error estima- 
tor shown in (9), which is based on the recovery of an enhanced stress field cr*. The 
accuracy of such estimates depends on the quality of the recovered field. To improve 
the quality of the recovered fields, numerical results indicate that the when recov- 
ering singular fields, statical admissibility and suitably chosen enrichment functions 
improve the accuracy [18, 35] 2 . In this work we consider the SPR-CX recovery tech- 
nique, which is an enhancement of the error estimator introduced in [26], to recover 
the solutions for both the primal and dual problems. The technique incorporates the 
ideas in [12] to guarantee locally on patches the exact satisfaction of the equilibrium 
equations, and the extension in [15] to singular problems. 

In the SPR-CX technique, as in the original SPR technique, a patch is defined 
as the set of elements connected to a vertex node. On each patch, a polynomial 
expansion for each one of the components of the recovered stress field is expressed 
in the form: 

p(x)a fc k = xx,yy,xy (45) 



x 



where p represents a polynomial basis and a are unknown coefficients. Usually, the 
polynomial basis is chosen equal to the finite element basis for the displacements. 

For the 2D case, the linear system of equations to evaluate the recovered stress 
field coupling the three stress components reads: 



<T*( X ) 






'p(x) 




a xx 1 


P(x)A = 


p(x) 


1 


a yy i 




p(x) 







(46) 



Constraint equations are introduced via Lagrange multipliers into the linear sys- 
tem used to solve for the coefficients A on each patch in order to enforce the satis- 
faction of the: 

• Internal equilibrium equations. 

• Boundary equilibrium equations: A point collocation approach is used to im- 
pose the satisfaction of a polynomial approximation to the tractions along the 
Neumann boundary intersecting the patch. 

• Compatibility equations: This additional constraint is also imposed to further 
increase the accuracy of the recovered stress field. 

Different techniques have been proposed to account for the singular part during 
the recovery process [15, 13, 16]. Here, following the ideas in [15], for singular 



2 The use of enrichment to improve recovery based error estimates for enriched approximations 
was discussed in detail in some of the first papers discussing derivative recovery techniques for 
enriched finite element approximations, see References [13, 14, 16]. A very detailed and clear 
discussion of a wide variety of error estimators and adaptive procedures for discontinuous failure 
is provided in [36]. 
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problems the exact stress field cr is decomposed into two stress fields, a smooth field 
<r smo and a singular field cr sing : 

C = CTsmo + CT sing- (47) 

Then, the recovered field cr* required to compute the error estimate given in (9) 
can be expressed as the contribution of two recovered stress fields, one smooth cr* mo 
and one singular cr* ing : 

Cr* = + cr s * ing . (48) 

For the recovery of the singular part, the expressions which describe the asymp- 
totic fields near the crack tip are used. To evaluate <x* ing we first obtain estimated 
values of the stress intensity factors K\ and K\\ using a domain integral method 
based on extraction functions [33, 37]. Notice that the recovered part cr* ing is an 
equilibrated field as it satisfies the equilibrium equations. 

Once the field cr* ing has been evaluated, an FE-type (discontinuous) approxima- 
tion to the smooth part cr^ mo can be obtained subtracting cr* ing from the raw FE 
field: 

ff L = ^ - ^sing- (49) 

Then, the field cr* mo is evaluated applying the enhancements of the SPR tech- 
nique previously described, i.e. satisfaction of equilibrium and compatibility equa- 
tions in each patch. Note that as both cr* mo and cr* ing satisfy the equilibrium equa- 
tions, cr* also satisfies equilibrium on each patch. A similar approach is followed in 
the case where an initial stress or strain is applied. We subtract the a priori known 
information (cr or e ) from the FE stress field, then we perform the recovery and 
finally add again cr or e . 

To obtain a continuous field, the recovered stresses cr* are directly evaluated at an 
integration point x through a partition of unity procedure [29], properly weighting 
the stress interpolation polynomials obtained from the different patches formed at 
the vertex nodes of the element containing x: 



n 



Cr*( X )=^iV J (x)cr}(x J ), (50) 
J=l 

where Nj are the shape functions associated with the vertex nodes n v . However, 
this process introduces a lack of equilibrium s = YTj=i ViVj • cr} when evaluating 
the divergence of the internal equilibrium equation, as explained in [26, 27]. 

This recovery technique, used to obtain the recovered field for the primal problem 
cr* is also used to evaluate cr*, considering different loading conditions. Two remarks 
are worth being made here. First, the analytical expressions that define the tractions 
t and body forces b for the dual problem are obtained from the interpretation 
of the functional Q(u) in terms of t d and h d , as seen in Section 4. Second, to 
enforce equilibrium conditions along the boundary of the domain of interest (Dol), 
we consider a different polynomial expansion on each side of the boundary and 
enforce statical admissibility of the normal and tangential stresses. Suppose that we 
have a patch intersected by r, such that Q e = Qi >e U Q 2 ,e for intersected elements, 
as indicated in Figure 1. To enforce equilibrium conditions along Tj we define the 
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stresses Jpj, o"h 2 \^i 011 eacn s ^ e °f the internal boundary. Then, the boundary 
equilibrium along Tj is: 



QCofiita - °"f2 2 |rj = t r . 



(51) 



where tr, = [t x t y } T are the prescribed tractions, e.g. those obtained from (27), and 
Q is the matrix form of Cauchy's law considering the unit normal n to Tj such that 



Q 



, Tly 



n 

n, 



(52) 





■Cv. 


Q 2 








< 


— F 








r, 



Figure 1: Equilibrium conditions along internal boundaries. 



6 Numerical results 

In this section numerical tests using 2D benchmark problems with exact solutions 
are used to investigate the quality of the proposed technique. The first problem 
has a smooth solution whilst the second is a singular problem. For both models we 
assume plane strain conditions. The /i-adaptive FE code used in the numerical 
examples is based on Cartesian meshes independent of the problem geometry [38] , 
with linear quadrilateral (Q4) elements 3 . To represent the domain geometry accu- 
rately, the integrals in elements cut by the boundary are restricted to the part of the 
element within the domain as in [40, 41]. The integration procedure in these ele- 
ments is based on the definition of triangular integration subdomains within Q e and 
aligning with the geometry of the domain. In elements cut by curved boundaries, 
these triangular subdomains can have curved boundaries. In that case, the actual 
geometry is reproduced using a transfinite mapping [42, 43]. Dirichlet boundary 
conditions are imposed using Lagrange multipliers following a procedure similar to 
that described in [44] and, more recently [39]. In cases where the geometry bound- 
ary cuts the mesh close to a node, preconditioning techniques are necessary, and the 
authors recommend the simple domain decomposition approach of [45] 4 . 

3 The interested reader is referred to the recent paper by Moumnassi and colleagues [39] which 
discusses recent advances in "ambient space finite elements" and proposes hybrid level set/FEMs 
able to handle sharp corners and edges. 

4 The interested reader may want to refer to the most insightful papers of [46] and [47, 48], the 
latter two being based on algebraic multi grid approaches. 
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To assess the performance of the proposed technique we consider the effectivity 
index of the error estimator 9 defined as the quotient of estimated error Q(e es ) in 
the quantity of interest over the exact error Q(e): 

Q(e es ) 

9 = W (53) 

We can also represent the effectivity in the Qol 9q i, defined as the corrected 
value of the Qol Q(u h ) using the error estimate, divided by the exact value Q(u): 

9qo1 - — W) — ' ( } 

The relative error in the Qol for the exact and estimated error are 

Q( ) _ JQ(e)| Q( v \Q(e es )\ ( , 

The distribution of the local effectivity index D is analysed at element level, 
following the definitions described in [12] for the error in energy norm, adapted here 
to the case of the error in Qol: 



D = 9 e — 1 if 9 e > 1 n( „e 

D = i-i if f~<i with B ' J ik^ (56) 

9 e v ' 

where superscript e denotes evaluation at element level. To evaluate Q(e e ) and 
Q{e e es ) we use (13, 17) locally in each element. Nonetheless, we should remark that 
this is only possible for some problems with analytical solutions as the exact value 
of a is unknown in the vast majority of cases, especially for the dual problem. 

Once error in the Qol has been estimated, the local error estimates in each 
element can be used to perform /i-adaptive analyses using similar techniques to 
those already available for the error in the energy norm. The refinement of the mesh 
using the error estimate as the guiding parameter considers an stopping criterion 
that checks the value of the estimated error against a prescribed or desired error. If 
the estimated error is higher than the desired error then the mesh is refined. Several 
procedures to perform the refinement are available in the literature. To define the 
size of the elements in the new mesh we follow the adaptive process described in 
[49, 50, 51] which minimises the number of elements in the new mesh for a given 
accuracy level. This criterion is equivalent to the traditional approach of equally 
distributing the error in each element of the new mesh as shown in [52, 53]. These 
criteria provide the size of the elements in the new mesh as a function of i) the ratio 
of the estimated error in energy norm in the current mesh to the desired error in the 
new mesh, and it) the estimated error in energy norm in each element, which always 
take non negative values. They cannot be directly used in goal-oriented adaptivity 
because the local contributions to the global error in the Qol, evaluated in each 
element using (17), can take negative values. Thus, for our implementation using 
/i-adaptive routines developed for the error in energy norm we prepare as input the 
square root of the absolute value of the error in the Qol - following the relation 
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from the expressions in (8, 13) - and the ratio of the estimated error in the Qol in 
the current mesh to the desired error in the new mesh. 

The refinement technique provides a distribution of the required new element 
sizes. These sizes are specified in each element of the current mesh, which will be 
recursively split into 4 new elements until the sizes of the elements are smaller than 
the required size. A maximum difference of only one refinement level will be allowed 
between adjacent elements. In these cases C° continuity will be enforced by means 
of the use of multipoint constraints [54, 55]. 



6.1 Problem 1: Thick-wall cylinder subjected to internal 
pressure. 

The geometrical model for this problem and the initial mesh are represented in 
Figure 2. Due to symmetry, only 1/4 of the section is modelled. The domain of 
interest (Dol) Qi is denoted by a dark square whereas the contours of interest Tj 
and T are the internal and external surfaces of radius a and b. Young's modulus is 
E = 1000, Poisson's ratio is v = 0.3, a = 5, b = 20 and the load P = 1. 

The exact solution for the radial displacement assuming plane strain conditions 
is given by 

Pil + v) ( „ „ x b 2 ' 



u ' {r) = E^vA r{l - 2v)+ ^) (57) 

where c = b/a, r = \J x 2 + y 2 and (j) = arctan(y/x). Stresses in cylindrical coordi- 
nates are 

P ( b 2 \ P ( b 2 \ P 

^) = ^3l(l-^j ^) = ^3l( 1 + ^) * z (r,9) = 2v— (58) 




Figure 2: Thick-wall cylinder subjected to an internal pressure. Model and analyt- 
ical solution, the domain of interest fL is indicated in dark. 



Several linear quantities of interest have been considered for this problem: the 
mean radial displacements along r o , the mean displacements u x in the domain of 
interest O i; the mean stresses a x in f2j and the mean value of the normal tractions 
along Tj. For the last case, instead of applying pressure P, we impose displacements 
u r (a) on this surface, see (57). 
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6.1.1 Problem l.a.: Mean radial displacements along r o 

Let Q(u) be the functional that evaluates the mean normal displacement u n along 
T such that: 

Q(u) =u n = ^- [ (Ru) T c u dr (59) 

where R is the standard rotation matrix for the displacements that aligns the coor- 
dinate system with the boundary T Q and c u = {1,0} T is the extraction vector that 
selects the normal component. Thus, the exact value of the Qol given by (57) for 
r = b is u n = 2.426 • KT 3 . 

Figure 3 shows the set of meshes resulting from the /i-adaptive process guided 
by the error estimate in this Qol. 




a) Mesh 1 b)Mesh 2 




c) Mesh 3 d)Mesh 4 

Figure 3: Problem l.a. Sequence of meshes for the mean displacement in a domain 
of interest r o . 

In Table 1 we show the results for the error estimate £spr-cx = £i from (18), 
evaluated using the SPR-CX recovery technique both in the primal and the dual 
problems, and the exact error Q(e). The recovery technique accurately captures the 
exact error and provides good effectivities, with values of 9 = 0.9473 for the coarsest 
mesh with 180 degrees of freedom (dof). The evolution of the effect ivity index for 
the proposed technique and the standard SPR is represented in Figure 4. In this 
case, SPR-CX converges faster to the optimal value 9 = 1 and shows a high level of 
accuracy (9 = 0.9928 for 3294 dof). 
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Table 1: Problem l.a. The table shows the values for the error estimate £spr-cx 
and effectivities. 



dof 


Q(e es ) = 


^SPR-CX 


Q(e) 




9 




180 


5.909651 


• 10~ 5 


6.238291 


10" 5 


0.94731888 


0.99864571 


322 


1.739913 


■ 10~ 5 


1.755027 


10~ 5 


0.99138817 


0.99993772 


937 


4.397256 


■ 10~ 6 


4.459314 


10~ 6 


0.98608343 


0.99997443 


3,294 


1.077612 


■ 10~ 6 


1.085411 


10~ 6 


0.99281470 


0.99999679 



9 



1 






0.95 






0.9 






0.85 




^-^SPR-CX 
^-<?SPR 



10 2 10 3 
dof 



Figure 4: Problem l.a. Evolution of the effect ivity index 9 considering locally 
equilibrated , £spr-cx, and non-equilibrated recovery, £spr- 

Note that in this case, the dual problem corresponds to a traction t = R T c/|r o | 
that represents a constant traction normal to the external boundary. Therefore, the 
solution to the dual problem can be evaluated exactly for this quantity of interest 
using the analytical solution for a cylinder under external pressure: 

where P Q represents the applied external pressure. As the exact solution for the 
dual problem is available, it is possible to evaluate the local effectivity D in (56) 
at element level to evaluate the local quality of the error estimate in the quantity 
of interest. Figure 5 shows the evolution of the mean absolute value m(|D|) and 
standard deviation cr(D) of the local effectivity. The ideal scenario is that both 
parameters are as small as possible and go to zero as we increase the number of 
DOFs. In the figure we see that the SPR-CX gives a better local estimation, with 
values of cr(D) and m(\D\) that are smaller than those for the SPR - for the mesh 
with approx. 3300 dof a(D) = 0.11 for the SPR-CX, compared to a(D) = 0.36 for 
the SPR. 

As the /i-adaptive procedures use local information to refine the mesh, providing 
an accurate local error estimate to the adaptive algorithm is of critical importance. 
For this reason, the local performance of the proposed technique indicates that 
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m(\D\) a(D) 




10 2 10 3 10 2 10 3 

dof dof 



Figure 5: Problem l.a. Evolution of the mean absolute value m(|D|) and standard 
deviation cr(D) of the local effectivity considering locally equilibrated , £spr-cx, 
and non-equilibrated recovery, £spr- 



the error estimator based on equilibrated recovered fields for the primal and dual 
problems is superior to the standard recovery techniques to guide the /i-adaptive 
process, even in cases in which the global effectivity is similar to the effectivity of 
the standard SPR. 

6.1.2 Problem l.b.: Mean normal traction t n along Ti 

Consider as quantity of interest the mean normal traction along the inner boundary 
of the cylinder in Figure 2. As previously mentioned, we consider an equivalent 
model where displacements are imposed along the inner boundary. The radial dis- 
placements to impose are evaluated from (57) such that they correspond to an 
internal pressure of P = 1. In this case, given the displacements we are interested 
in evaluating the tractions as Qol. 

To formulate the dual problem we consider displacement constraints along Tj 
such that 

u r {r = a) = — (61) 

I i\ 

Table 2 shows the results for the error estimate £spr_cx using the proposed 
recovery technique and the exact error Q(e). The procedure accurately captures 
the error and yields good effectivities for this Qol, with values of 9 = 0.9109 for the 
last mesh with 2218 dof. Figure 6 shows the evolution of the effectivity with the 
number of dof for the SPR-CX and the SPR. Results show a better performance of 
the proposed technique when compared with the SPR which presents an oscillatory 
behaviour (9 = 0.9466 for mesh 2 with 260 dof and 9 = 0.5681 for mesh 3 with 732 
dof). 

Figure 7 represents the evolution of the mean absolute value m(\D\) and standard 
deviation a(D) of the local effectivity. Again, for this example, the SPR-CX gives 
lower values of these parameters than the SPR. The improved local performance of 
the SPR-CX is particularly useful for the adaptive algorithm. Notice that for the 
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Table 2: Problem l.b. The table shows the values for the error estimate £spr-cx 
and effectivities. 



dof 


Q(e es ) = 


^SPR-CX 


-Q(e 


) 




9 


&QoI 


164 


2.245490 


• io- 1 


4.242762 • 


10" 


-1 


0.52925199 


1.64734771 


256 


1.794064 


■ io- 2 


2.078322 • 


10" 


-2 


0.86322691 


1.03855364 


652 


3.910983 


• IO" 3 


4.170597- 


10" 


-3 


0.93775134 


1.00808981 


2,218 


1.008523 


• IO" 3 


1.107166- 


io- 


-3 


0.91090543 


1.00211542 



9 





-H 1 1 1 — 




-©-^SPR-CX " 




-■^^SPR 



10 2 10 3 
dof 



Figure 6: Problem l.b. Evolution of the effectivity index 9 considering locally 
equilibrated , £spr-cx, and non-equilibrated recovery, £spr- 

second mesh, although the global effectivity is close to unity for the SPR, m{\D\) 
and cr(D) indicate that the SPR-CX is a superior choice when estimating the true 
error at the element level. Thus, the apparent satisfactory behaviour of the SPR in 
mesh 2 is due to error compensations of large values of local effectivities. 

6.1.3 Problem I.e.: Mean displacements u x in fl f 

Let us consider the mean displacement u x in Qj as the quantity of interest. The 
objective is to evaluate the error when evaluating u x defined by the functional: 

Q(u) = u x = ^- I u x dn (62) 

The exact value of the Qol can be computed for this problem and is u x = 
0.002238239291713. Figure 8 shows the first four meshes used in the /i-adaptive 
refinement process guided by the error estimate for the Qol. 

Figure 9 shows the relative error (in percentage) for the error estimates in (18) 
evaluated using the proposed recovery technique. The exact error Q(e) is also repre- 
sented for comparison. For all the curves the relative error decreases monotonically 
when increasing the number of dof, indicating that the /i-adaptive process has a 
stable convergence. The most accurate estimation is obtained for the estimate Si, 
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m(\D\) 



a(D) 




10- 1 



10 1 



10 U : 





^^^SPR-CX 




^-^SPR 





i ii 



10 2 



10 3 



dof 



Figure 7: Problem l.b. Evolution of the mean absolute value m(|Z)|) and standard 
deviation &(D) of the local effectivity considering locally equilibrated , £spr-cx, 
and non-equilibrated recovery, £spr- 
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a) Mesh 1 



b)Mesh 2 





c) Mesh 3 d)Mesh 4 

Figure 8: Problem I.e. Sequence of /i-adaptive refined meshes. 



which practically coincides with the exact relative error. The other estimates tend to 
overestimate the exact error, although strictly speaking they do not have bounding 
properties, which is in agreement with the observations made previously in Sec- 
tion 3.2. Figure 10 represents the evolution of the effectivity index as we increase 
the number of dof. The curves are in consonance with the results in Figure 9 and 
show good values for the effectivity index, especially for E\ with values between 
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[0.9898,1.0481]. 

T) Q {%) 




10 3 10 4 10 5 

dof 



Figure 9: Problem I.e. Evolution of the relative error rp considering the error 
estimates in (18) and the exact error Q(e). 




Figure 10: Problem I.e. Evolution of the effectivity index 9 for the error estimates 
in (18) 

Table 3 shows the estimated error in the Qol Q(e es ), the exact error Q(e), the 
effectivity in the quantity of interest 9q j and the effectivity of the error estimator 
9 for the sharp error estimate £spr-cx = £i- Comparing Q(e) and Q{e es ) we can 
notice that both values decrease as we refine and that the estimate £spr-cx gives 
a good approximation of the exact error. The effectivity of the error estimator 9 
converges and is very close to the theoretical value 9 = 1 (with 9 = 1.0123 for 75323 
dof). As expected from these results, the effectivity 9q q i is very accurate as well, 
with 9 QoI = 1.00000012 for 75323 dof. 

If in (18a) we consider the case of the non-equilibrated super convergent patch 
recovery procedure, resembling the averaging error estimators presented in [24], we 
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Table 3: Problem I.e. Error estimate £spr-cx and its effectivities. 



dof 


Qi^es) — £sPR- 


-cx 


Q(e) 






9 


9qoi 


528 


1.670429 • 


10" 


-5 


1.666535 ■ 


io- 


-5 


1.00233710 


1.00001740 


1,126 


1.849599 ■ 


10" 


-6 


1.868635 • 


io- 


-6 


0.98981317 


0.99999150 


4,348 


4.240859 • 


10" 


-7 


4.046075 • 


io- 


-7 


1.04814162 


1.00000870 


18,076 


1.017070 • 


10" 


-7 


9.922496 • 


io- 


-8 


1.02501376 


1.00000111 


75,328 


2.128196 • 


io- 


-8 


2.102331 • 


io- 


-8 


1.01230285 


1.00000012 



obtain the results shown in Figure 11. This figure shows that, in this case, the 
effectivity of the error estimator provided by the standard SPR technique is similar 
to the effectivity obtained with the SPR-CX technique here proposed, although the 
latter results in more accurate values for the coarsest mesh (9 = 1.1659 for the SPR 
and 9 = 1.0023 for the SPR-CX, considering 528 dof). In any case, we should recall 
that the local behaviour with the SPR-CX tends to be better than with the SPR. 

9 



1.2 



1 



10 3 10 4 10 5 

dof 

Figure 11: Problem I.e. Evolution of the effectivity index 9 considering equilibrated 
^spr-cx and non-equilibrated recovery, £spr. 



6.1.4 Problem l.d.: Mean stress a x in f2, 

Consider as Qol the mean stress value a x given in (31). Figure 12 shows the first 
four meshes of bilinear elements used in the refinement process guided by the error 
estimated for this Qol. 

The evolution of the relative error for the estimates presented in (18) using the 
proposed recovery technique and the exact error is shown in Figure 13. Similar to the 
observations done for the previous examples, the most accurate results are obtained 
when considering the estimate E\. In this case, the other estimates considerably 
overestimate the true error. Figure 14 shows the effectivity index for £spr-cx, 
which uses the locally equilibrated SPR-CX recovery technique, together with the 
effectivity obtained with the non-equilibrated SPR technique (£spr curve). This 




23 




graph clearly shows the improvement obtained with the SPR-CX recovery technique, 
with effectivities very close to 1 (6 = 0.9797 for the SPR-CX whilst 9 = 1.6209 for 
the SPR, considering 19573 dof), in contrast with the unstable behaviour of the 
SPR technique. 




10 3 10 4 10 5 

dof 



Figure 13: Problem l.d. Evolution of the relative error rp considering the error 
estimates in (18) and the exact error Q(e). 
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1 







10 3 10 4 10 5 

dof 

Figure 14: Problem l.d. Evolution of the effectivity index 9 considering equilibrated, 
£spr-cx, and non-equilibrated recovery, £spr. 

Table 4 shows the estimate Q(e es ) = £spr-cx, the exact error Q(e), the effec- 
tivity index for the Qol 9q q j and the effectivity of the error estimator 9. For this 
problem the exact value of the Qol is a x = 0.06. Table 4 indicates that the equili- 
brated recovery procedure (SPR-CX) provides very accurate estimations of the error 
in the Qol and in the value of the Qol itself, with 9 = 1.0259 and 9 QoI = 1.00000063 
for 79442 dof. 

Table 4: Problem l.d. Values for the error estimate £spr_cx and effectivities. 



dof 




^SPR-CX 


Q(e) 






9 


®QoI 


528 


1.659077 


• 10" 


-3 


1.738923 


io- 


-3 


0.95408307 


1.00119769 


1,350 


8.225557 


• 10" 


-5 


7.687113 


io- 


-5 


1.07004501 


1.00008077 


5,131 


2.678276 


• io- 


-5 


3.016986 


io- 


-5 


0.88773219 


0.99994919 


19,573 


1.108891 


• io- 


-5 


1.131773 


io- 


-5 


0.97978268 


0.99999657 


79,442 


1.655151 


• io- 


-6 


1.613254 


io- 


-6 


1.02597049 


1.00000063 




6.2 Problem 2: L-Shape plate 

Let us consider the singular problem of a finite portion of an infinite domain with 
a reentrant corner. The model is loaded on the boundary with the tractions 
corresponding to the first terms of the asymptotic expansion that describes the 
exact solution under mixed mode loading conditions around the singular vertex, see 
Figure 15. The exact values of boundary tractions on the boundaries represented 
by discontinuous thick lines have been imposed in the FE analyses. 

The exact displacement and stress fields for this singular elasticity problem can 
be found in [33]. Exact values of the generalised stress intensity factors (GSIF) 
[33] under mixed mode have been taken as Kj — 1 and Ku = 1. The material 
parameters are Young's modulus E = 1000, and Poisson's ratio v = 0.3. As 
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Figure 15: Problem 2. L-shaped domain. 



the analytical solution of this problem is singular at the reentrant corner of the 
plate, for the recovery of the dual and primal fields we apply the singular + smooth 
decomposition of the stresses as explained in Section 5. We utilise a domain integral 
method based on extraction functions to obtain an approximation of the recovered 
singular part as explained in Section 5. 

In this example, we consider the GSIFs K\ and K\\ as the quantities of interest. 
Figure 16 shows the Cartesian meshes used to solve the primal and dual problems. 
For the dual problem, we use the same Dirichlet conditions as shown in Figure 15 
and the set of nodal forces used to extract the Qol in the annular domain defined 
by a plateau function q, shown in Figure 16. Function q is defined such that q = 1 
for r < 0.6, q = for r > 0.8 and has a smooth transition for 0.6 < r < 0.8 given 
by a quartic spline 1 — 6s 2 + 8s 3 — 3s 4 . 




Figure 16: Problem 2. Cartesian meshes with /i-adaptive refinement. 



In order to impose equilibrium conditions during the recovery of the stresses 
we use the following approach. For the primal solution, in each patch, we enforce 
internal equilibrium and compatibility in Q, and boundary equilibrium all along the 
Neumann boundary. For the dual problem, we enforce internal equilibrium using the 
initial strains and body loads shown in (44) and homogeneous Neumann boundary 
conditions. Compatibility is enforced within all the domain. 



26 



Table 5: Problem 2. Stress intensity factor K\ as Qol using initial strains formula- 
tion. 



dof Q(e es ) = £spr-cx Q(e) 9 9 



Qol 



1,101 7.302834 • 10~ 3 9.566751 • 10~ 3 0.76335570 0.99773608 

2,927 2.065838 • 10~ 3 2.206343 • 10" 3 0.93631774 0.99985950 

10,399 5.194502 • 10" 4 5.280162 • 10" 4 0.98377707 0.99999143 

39,193 1.300133 • 10" 4 1.307862 • 10~ 4 0.99409058 0.99999923 



Table 6: Problem 2. Stress intensity factor Ku as Qol using initial strains formula- 
tion. 



dof 


Q(e es ) = 


^SPR-CX 


Q(e) 






9 


9qoi 


1,101 


2.031449 


• io- 3 


1.797716 • 


10" 


-3 


1.13001659 


1.00023373 


2,927 


5.273469 


• io- 4 


4.342458 • 


10" 


-4 


1.21439713 


1.00009310 


10,399 


1.013432 


• 10~ 4 


9.947112 • 


10" 


-5 


1.01881993 


1.00000187 


39,193 


2.440128 


• 10~ 5 


2.387141 • 


io- 


-5 


1.02219675 


1.00000053 



Table 5 shows the results for the stress intensity factor K\. Similarly to the 
results for other Qols, we observe that the proposed technique provides an accurate 
estimate Q(e es ) of the exact error Q(e). The effectivity index 9 is always close to 
the optimal value 9 = 1 and for the last mesh with 39,193 dof we obtain 9 = 0.9940. 
As a result, and in agreement with the previous cases, the effectivity in the Qol is 
highly accurate, with 9q q i = 0.9999 for the same mesh. Table 6 shows the same 
results for the stress intensity factor Ku- Again, we observe a satisfactory behaviour 
of the error indicator and very accurate effectivities, both for the error estimate and 
for the Qol itself (9 = 1.0222 and 9 QoI = 1.0 for the last mesh). 

Figure 17 shows the evolution of the relative error with respect to the number of 
dof for the two Qols. Figures 18 and 19 show the evolution of the effectivity index of 
the error estimators obtained with the locally equilibrated SPR-CX technique, and 
with the non-equilibrated recovery as we increase the number of dof. The results 
indicate that the proposed methodology accurately evaluates the error in the Qol, 
giving values of 9 close to 1 and considerably improving the results obtained with 
the original SPR technique (for the finest mesh considering Ku we have 9 = 0.8610 
for the standard SPR and 9 = 1.0222 for the SPR-CX). 



7 Conclusions and future work 

In this paper, we present an a posteriori recovery-based strategy that aims to control 
the error in quantities of interest. The proposed technique is based on the use of 
a recovery procedure that considers locally equilibrated stress fields for both the 
primal and the dual problem. 

To recover the solution for the primal problem the formulation considers the 
enforcement of equilibrium and compatibility equations, and the splitting of the FE 
field for singular problems [15]. In order to enforce equilibrium conditions when 
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Figure 17: Problem 2. Evolution of the relative error 77^ for Qols Kj and Ku, for 
the exact and estimated error as defined in (55). 
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Figure 18: Problem 2. Evolution of the effect ivity index 9 for K\ 

evaluating the recovered stress fields for the dual problem, we have considered an 
approach that allows to express the functional which defines some Qols in terms 
of body loads and boundary tractions applied to the dual problem. The proposed 
technique was tested on different quantities of interest: mean displacements and 
stresses on a domain of interest, mean displacements and tractions along a boundary 
and the generalised stress intensity factor in the singular problem of a reentrant 
corner. 

The methodology we proposed provides accurate global and local evaluations 
of the error in the different quantities of interest analysed, improving the results 
obtained with the original SPR technique, whose satisfactory estimation of the global 
error can be attributed to compensating errors. In particular, our approach proves 
superior to the standard SPR when considering the local error. 

An extension of the work presented here for extended finite element approxima- 
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Figure 19: Problem 2. Evolution of the effectivity index 9 for K\\ 

tions is now under development. The final goal is to be able to guide the adaptive 
process in 3D XFEM problems modelling fatigue crack growth using the stress in- 
tensity factors as design parameters for industrial applications of the XFEM similar 
to those discussed in the early work of References [56] [57], [58] and [59, 60]. 
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